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The question of whether glass continues to relax at low temperature is of fundamental and practical 
interest. Here, we report a novel atomistic simulation method allowing us to directly access the long¬ 
term dynamics of glass relaxation at room temperature. We find that the potential energy relaxation 
follows a stretched exponential decay, with a stretching exponent /3 = 3/5, as predicted by Phillips’ 
diffusion-trap model. Interestingly, volume relaxation is also found. However, it is not correlated to 
the energy relaxation, but is rather a manifestation of the mixed alkali effect. 

PACS numbers: 65.60.+a, 62.40.+i, 81.05.Kf 


While it is indeed commonly believed that, as frozen 
supercooled liquids, glasses should continue to flow over 
the years (e.g. in the case of the stained-glass windows of 
medieval cathedrals mm ), the dramatic increase of their 
viscosity below the glass transition temperature T g sug¬ 
gests, on the contrary, that their relaxation time is on the 
order of 10 32 years at room temperature [3]. However, a 
recent study conducted by Mauro et al. H] reported the 
intriguing dynamics of the relaxation of a commercial 
Corning® Gorilla® Glass at room temperature, over 1.5 
years. Fused silica [5] and metallic glasses [6] were also re¬ 
ported to show an appreciable relaxation over time. More 
generally, the relaxation of glasses at low temperature is 
known as the "thermometer effect", originating from the 
fact that, in the 19 st century, thermometers made of al¬ 
kali lime silicate glass were experiencing gradual changes 
of dimension over time, making them inaccurate IB US¬ 
As it cannot be described as a viscous process, the 
physical origin of such room temperature relaxation re¬ 
mains unclear. However, the shape of the relaxation 
function contains valuable information to discriminate 
between the different proposed models. In particular, 
Mauro et al. observed that the volume of the relax¬ 
ing Gorilla® Glass followed a stretched exponential de¬ 
cay function [ 3 ], with a stretching exponent /3 = 3/7, 
which, interestingly, is the value predicted by the Phillips 
diffusion-trap model for a relaxation dominated by long- 
range pathways [ 5 ]. Understanding and predicting the 
relaxation dynamics of glasses has important practical 
applications, e.g. for optical fibers m, substrate glass 
for liquid crystal displays m , and chemically strength¬ 
ened cover glass for smartphones and tablets m- 

Here, based on molecular dynamics (MD) simulations, 
we investigate the atomistic origin of the relaxation of 
realistic alkali aluminosilicate glasses at temperature far 
below T g . As traditional MD is unable to reproduce year¬ 


long relaxations, we introduce a novel artificial relaxation 
method based on cyclic stress perturbations. The poten¬ 
tial energy of the simulated glasses is found to follow a 
stretched-exponential decay function, with a stretching 
exponent /3 = 3/5, characteristic of systems dominated 
by short-range forces |9j. On the contrary, the volume 
relaxation appears composition-specific and strongly cor¬ 
related with the dynamics of the alkali atoms. 

For this study, we simulated a 2991 atoms mixed potas¬ 
sium sodium aluminosilicate glass (denoted KNAS here¬ 
after), of composition (K 2 0) 8 (Na20) 8 (Al20 3 )9(Si02)75, 
which we expect to be representative of that of the per- 
alkaline Corning® Gorilla® Glass used in Ref. |!J. In 
addition, in order to evaluate the effect of the mixed 
alkali effect on relaxation, we simulated a similar sin¬ 
gle alkali glass (denoted NAS hereafter) of composition 
(Na20)i6(Al 2 03)9(Si0 2 )75. All MD simulations are per¬ 
formed with the LAMMPS package im using the well- 
established Teter potential EHH! and an integration 
time-step of 1 fs. Coulomb interactions were evaluated by 
the Ewald summation method, with a cutoff of 12 A. The 
short-range interaction cutoff was chosen at 8.0 A. Liq¬ 
uids were first generated by placing the atoms randomly 
in the simulation box. The liquids were then equilibrated 
at 5000 K in the NPT ensemble (constant pressure) for 1 
ns, at zero pressure, to assure the loss of the memory of 
the initial configuration. Glasses were formed by linear 
cooling of the liquids from 5000 to 300 K with a cooling 
rate of 1 K/ps in the NPT ensemble at zero pressure. 
The obtained glasses were eventually relaxed for an ad¬ 
ditional 1 ns in the NPT ensemble at 300 K and zero 
pressure. A detailed structural analysis of the simulated 
glasses can be found in Ref. El- 

Traditional MD simulations are usually limited to a 
few nanoseconds, which prevents one from using them to 
predict long-term relaxation at low temperature. On the 
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other hand, kinetic Monte-Carlo simulations like ART 
m offer an attractive alternative to perform simulations 
up to a few seconds, but their application to alkali silicate 
is challenging due to the high mobility of the alkali atoms, 
which results in a huge number of small energy barriers 
to compute. Since a direct simulation of relaxation dy¬ 
namics is, at this point, unachievable, we developed a 
new method to simulate an artificial relaxation. 

This method is inspired by the artificial aging or reju¬ 
venation observed in granular materials subjected to vi¬ 
brations DM- Small vibrations induce a compaction 
of the material, that is, they make the system artificially 
age. On the other hand, large vibrations randomize the 
grain arrangements, which decreases the overall compact¬ 
ness and, therefore, make the system rejuvenate. Similar 
ideas, relying on the energy landscape approach mma, 
have been applied to amorphous solids, based on the fact 
that small stresses deform the energy landscape under¬ 
gone by the atoms. This can result in the removal of 
some energy barriers existing at zero stress, thus allow¬ 
ing atoms to jump over them in order to relax to lower 
energy states (see the schematic representation in Fig. 
[TJ. This transformation is irreversible as, once the stress 
is removed, the system remains in its "aged" state. On 
the contrary, large stresses move the system far from its 
initial state, which eventually leads to rejuvenation, sim¬ 
ilar to thermal annealing (2?, S3]. As such, several ex¬ 
perimental and simulation studies reported that glasses 
can undergo overaging or rejuvenation, when subjected 
to small or large shear stresses, respectively |21fl2i?] . The 
novel method presented here is based on a succession of 
many of such stress perturbations, inducing successive 
small relaxations of the system. 

In practice, a glass initially relaxed at zero pressure is 
subjected to athermal stress perturbations cycles. Each 
perturbation cycle consists, successively, of (1) an hydro¬ 
static compressive stress —<to and (2) an hydrostatic ten¬ 
sile stress (To- Therefore, the average stress remains zero 
over each cycle. During both of the phases (1) and (2), 
a volume-variant energy minimization is performed, in 
which the positions of the atoms and the volume/shape 
of the simulation box are updated in an iterative pro¬ 
cess, in order for the system to achieve simultaneously 
both an energy minimum for the potential energy of 
the atoms U and a stress er close to ±<7o- This com¬ 
bined optimization is performed by minimizing a ficti¬ 
tious energy E^ c = U + E s tra in, where U is the potential 
energy and if s train a strain energy expression proposed 
by Parrinello and Rahman [30]. Hence, during each cy¬ 
cle, the system is allowed to jump over energy barriers, 
with a roughly constant activation energy on the order of 
-Estrain- According to transition state theory EL each cy¬ 
cle should roughly correspond to a constant time of about 
At = Vq 1 exp(-Estram/& bT), where u 0 is a characteristic 
frequency of attempt, ks the Boltzmann constant, and T 
the temperature [32]. Therefore, the relaxation observed 



FIG. 1. (color online). Relative variation of the poten¬ 
tial energy (stabilization) of the simulated mixed potassium 
sodium aluminosilicate (KNAS) and sodium aluminosilicate 
(NAS) glasses, with respect to the number of stress perturba¬ 
tion cycles applied N (with an stress amplitude of 0.4 GPa). 
For better clarity, the curve obtained for the NAS glass is 
shifted by —1. The simulated data are fitted by stretched 
exponential decay functions with different stretching expo¬ 
nents (3 (solid lines). The inset shows the relative variation 
of the potential energy of the simulated KNAS glass, for dif¬ 
ferent amplitudes of stress perturbations. The cartoon on the 
bottom-left is a schematic representartion of a strain-induced 
disappearance of energy barrier. The W curves represent the 
potential energy between two meta-stable equilibrium states, 
and the circles represent the state of the system. The arrow 
corresponds to a strain-activated jump of the system towards 
the more stable equilibrium state. 


with respect to the number of stress perturbations cy¬ 
cles applied N should be representative of the relaxation 
versus time t, with a converting factor t = N At. 

Figure [I] shows the relative variation of the potential 
energies of the KNAS and NAS glasses, with respect to 
the number N of stress perturbation cycles applied, with 
a sub-yield amplitude cto = 0.4 GPa. As expected, the 
stress perturbations allow the system to relax towards 
lower energy states. This stabilization is gradual and 
about 10 5 of such cycles are needed for the potential 
energy to plateau. Eventually, a significant decrease of 
energy (around 0.4%) is achieved. We note that the en¬ 
ergy relaxation profile appears to be very similar for the 
KNAS and NAS glasses. The influence of <to was as¬ 
sessed by performing similar simulations with cr 0 = 0.7, 
1, 2, and 5 GPa. As shown in the inset of Fig. [T] no 
significant change of the shape of the energy relaxation 
is observed at low <7q. Therefore, the relaxation dynam¬ 
ics does not depend on the choice of <Jo- However, the 
energy starts to decrease more slowly for a 0 = 5 GPa, 
which suggests that, for such a high stress, the system is 
slightly rejuvenated at each cycle. In the following, we 
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FIG. 2. (color online). Linear strain (compression when neg¬ 
ative) observed during the low temperature relaxation of the 
simulated mixed potassium sodium aluminosilicate (KNAS) 
and sodium aluminosilicate (NAS) glasses (left axis) with re¬ 
spect to the number of stress perturbation cycles applied N. 
Simulated results are compared to the experimental measure¬ 
ment of the room temperature compaction of a commercial 
KNAS glass (right axis) with respect to the time of relax¬ 
ation [1]. Experimental data are fitted by a stretched expo¬ 
nential decay function (solid line) with a stretching exponent 
P = 3/7. 


keep do = 0.4 GPa. Overall, this result constitutes the 
first direct simulation of the relaxation of a glass at low 
temperature, to the best of our knowledge. 

Interestingly, the volume of the glasses does not remain 
constant through the relaxation procedure. As shown 
in Fig. [2] the simulated KNAS glass shows a gradual 
compaction with respect to the number of stress pertur¬ 
bation cycles applied, eventually achieving a final linear 
strain of about —0.5%. As the average applied stress 
remains zero, this compaction cannot be explained by 
elastic deformations. We note that the extent of the com¬ 
paction is higher in the present simulation than observed 
experimentally |1J. This discrepancy is likely to arise 
from the cooling rate used in simulation, which is much 
higher than the ones typically achieved experimentally, 
and from the different compositions of the glasses. How¬ 
ever, as shown in Fig. [2] the general trend of the volume 
relaxation is remarkably fairly similar to that observed 
experimentally. 

The origin of the densification of the Gorilla® Glass 
and, more generally, of the thermometer effect, has been 
attributed to the mixed alkali effect (MAE) [7] Hj. An 
MAE is generally observed in silicate or borate glasses 
containing at least two different network modifying alkali 
oxides AO 2 and BO 2 . It manifests by a strong non linear 
evolution of the transport properties with respect to the 
fraction A/B. Here, the compaction of the KNAS glass 
can be clearly attributed to the MAE. Indeed, as shown 



FIG. 3. (color online). Relaxation function / for the energy 
and volume, computed during the relaxation of the simulated 
mixed potassium sodium aluminosilicate (KNAS) glass, with 
respect to the normalized number of stress perturbation cy¬ 
cles applied N/No (see text). The computed data are fitted 
with stretched exponential decay functions, with stretching 
exponents /3 = 3/5 and 1 for the energy and the volume, 
respectively. The inset shows log(—log(/)) with respect to 
\og(N/No), for the energy and volume, whose slope is equal 
to /3. Solid lines with the slopes /3 = 3/7, 3/5, and 1 are 
added for comparison. 


in Fig. [2] the volume of the single alkali NAS glass fea¬ 
tures a completely different trend, remaining fairly con¬ 
stant over time. This result is consistent with the ob¬ 
servation that single alkali-containing Gorilla® Glass 2 
and 3 do not show any discernable volume relaxation [3] . 
These results highlight the fact that volume relaxation is 
not an intrinsic feature of glass, but strongly depends on 
composition. On the contrary, based on the present re¬ 
sults, energy relaxation appears to be more general and, 
therefore, more relevant to characterize the intrinsic dy¬ 
namics of glass relaxation. 

Indeed, the general shape of the relaxation curve of the 
potential energy is of fundamental importance as it man¬ 
ifests from the topology of the relaxation process [33]. 
Glass relaxation behaviors are generally found to follow 
a stretched exponential decay function /(f) [53] : 

m= ex p(-(t/r)f > ) (i) 

where r is the relaxation time, and /3 a dimensionless 
stretching exponent satisfying 0 < /3 < 1. The case /3 = 1 
corresponds to a simple exponential decay. The diffusion- 
trap model from Phillips [9] , based on diffusion of excita¬ 
tions to randomly distributed traps, predicts that a theo¬ 
retical value for the stretching exponent /3 = d*/(d* + 2), 
where d* is the effective dimensionality of the channels 
along which the excitations diffuse in the configuration 
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space. This is described as d* = fd , where d is the di¬ 
mensionality of the system (3 for structural glasses), and 
/ is the proportion of the active channels of relaxation. 
Hence, when all the channels are active (/ = 1), one gets 
P = 3/5. When only long-range channels are active, by 
assuming an equipartitioning of the short- and long-range 
contributions (/ = 1/2), the model predicts (3 = 3/7. 

As shown in Fig. [2] the experimental volume relax¬ 
ation of the Gorilla® Glass was found to follow such a 
stretched exponential decay function, with a relaxation 
time t = 27.6 days and a stretching exponent /3 = 3/7, 
attributed to a relaxation dominated by long-range path¬ 
ways. Interestingly, the energy relaxation predicted by 
the present simulations also features a stretched exponen¬ 
tial decay, with an equivalent relaxation time r = Ay At, 
with No = 3200 and 1300 for the KNAS and NAS glasses, 
respectively. This allows us to roughly evaluate At, the 
fictitious time elapsed during each stress perturbation cy¬ 
cle, to be on the order of 10 min. However, as shown in 
Fig. 0 we unambiguously find a stretching exponent 
fj = 3/5, which, as mentioned previously, corresponds 
to the situation in which all the relaxation channels are 
active [9] , which is typically observed for the stress relax¬ 
ation in glasses [33J [33] ■ We note that the /? = 3/5 factor 
offers a good fit, both for the KNAS and NAS glasses, 
which suggests that it does not depend on changes of 
composition. 

Surprisingly, as shown in Fig. [3j the energy and vol¬ 
ume relaxations of the KNAS glass do not show the same 
shape. First, we find Nq = 3200 and 12000 for the en¬ 
ergy and volume, respectively, suggesting a slower relax¬ 
ation for the volume than for the energy. Second, the 
volume relaxation actually follows a simple (that is, not 
stretched) relaxation decay function. As shown in the 
inset of Fig. [3] this difference of stretching exponent 
can be clearly established by plotting log(—log(/)) ver¬ 
sus log(A r /A r 0 ), whose slope is equal to /3. This result 
suggests that energy and volume relax in a fundamentally 
different way. More specifically, the relaxation of energy 
can be explained by the diffusion-trap model with all re¬ 
laxation channels being active [5] , whereas the relaxation 
of volume, showing a simple exponential relaxation, can 
be explained by a simple two-state model [36] . 

In order to elucidate the origin of the simple exponen¬ 
tial decay for the volume, we investigated the relaxation 
dynamics of each atomic species by calculating the inter¬ 
mediate scattering function (ISF) F s (k , N), computed at 
fc max , the position of the first peak of the structure factor 
m Interestingly, the ISF of Na and K atoms follow a 
trend similar to that of the volume relaxation. Indeed, as 
shown in Fig. [4] it can be fitted by a simple exponential 
decay function (/3 = 1) and show relaxation times similar 
to that of the volume relaxation (e.g., No = 11900 for Na 
atoms, as compared to 12000 for the volume relaxation). 
This clearly establishes that the volume relaxation ob¬ 
served in mixed alkali aluminosilicate glasses arises from 



FIG. 4. (color online). O, Si, Al, Na, and K intermediate 
scattering function F a (k,N) (computed at fc max , the position 
of the first peak of the structure factor) of the simulated potas¬ 
sium sodium aluminosilicate (KNAS) glass with respect to the 
number of stress perturbation cycles applied N. The simu¬ 
lated data are fitted by simple exponential decay functions 
(/3 = 1, solid lines). The inset shows the self-part of the van 
Hove correlation function 4nr 2 G B (r, N) for O, Si, Al, Na, and 
K, with respect to the distance r, computed at IV max = 2 16 
cycles. 


the relaxation dynamics of the alkali atoms. The self part 
of the van Hove correlation function, shown in the inset of 
Fig- ED shows that this relaxation of the alkali atoms con¬ 
sists of jumps, with a typical hopping distance of around 
3 A |35i - PTO] . On the contrary, the network former atoms 
show very limited motion, lower than 1 A. This highlights 
the fact that the volume relaxation is not correlated to 
the viscosity of the glass network, but rather to the dy¬ 
namics of the alkali modifiers only. The fact that, due to 
the MAE effect, this dynamics strongly depends on the 
composition of the glass may explain the discrepancy of 
stretching exponent between the simulation and experi¬ 
mental results. 

Featuring a stressed exponential decay function with 
/3 = 3/5, the shape of the energy relaxation is, accord¬ 
ing to the present results, more intrinsic to the glassy 
state than that of the volume relaxation. However, if 
the shape appears universal, the typical relaxation time 
r strongly depends on composition. In particular, relax¬ 
ation has been shown to be fastest for isostatic systems, 
which are characterized by an atomic network rigid but 
free of internal stress 1371 SI]. Beyond glasses, being able 
to predict and tune the relaxation and aging of materi¬ 
als could improve the understanding of memory encoded 
materials [52] [53] or protein folding [55155] . 
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